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ABSTRACT 

Using a set of high-resolution A/-body/SPH cosmological simulations with identical initial conditions but 
run with different numerical setups, we investigate the influence of baryonic matter on the mass distribution 
of dark halos when radiative cooling is not included. We compare the concentration parameters of about 400 
massive halos with virial mass from 10 13 to 7.1 x 10 14 h M s . We find that the concentration parameters for 
the total mass and dark matter distributions in nonradiative simulations are on average larger by ~ 3% and 
10% than those in a pure dark matter simulation. Our results indicate that the total mass density profile is little 
affected by a hot gas component in the simulations. After carefully excluding the effects of resolutions and 
spurious two-body heating between dark matter (DM) and gas particles we conclude that the increase of the 
DM concentration parameters is due to interactions between baryons and DM. We demonstrate this with the 
aid of idealized simulations of two-body mergers. The results of individual halos simulated with different mass 
resolutions show that the gas profiles of densities, temperature and entropy are subjects of mass resolution of 
SPH particles. In particular, we find that in the inner parts of halos, as the SPH resolution increases the gas 
density becomes higher but both the entropy and temperature decrease. 
Subject headings: galaxies: halos - cosmology: theory - dark matter -methods: numerical 



1. INTRODUCTION 

One potential problem of the cold dark matter (CDM) cos- 
mology is the mass distributions in the central region of clus- 
ters of galaxies. There is yet no consensus emerging from 
the observations. In general, simulations produce steep inner 
density profiles for DM halos (Navarro et al. 1997; Moore 
et al. 1998; Jing & Suto 2000), while observations of some 
clusters seem to prefer flat, core-like profiles (e.g., Tyson et 
al. 1998), and others prefer cusped profiles (e.g., Lewis et al. 
2003; Buote & Lewis 2004; Pointecouteau et al. 2005). Sand 
et al. (2004) show that the central DM profiles of clusters 
of galaxies have a slope of about -0.5, which is substantially 
flatter than the inner slope (-1) of the Navarro-Frenk- White 
(NFW) profile as found in CDM simulations. However, other 
authors argued that their interpretations are compromised by 
the assumption of sphericity (e.g., Bartelmann & Meneghetti 
2004; Meneghetti et al. 2005; Dalai & Keeton 2003). Re- 
cent observations of strong gravitational lensing have shown 
that some clusters of galaxies have very high mass concentra- 
tion parameters (e.g., Broadhurst et al. 2005) compared with 
the average value obtained in the DM simulations. The num- 
ber of giant arcs predicted by A/-body simulations may also 
be somewhat too low to be compatible with observations (Li 
etal. 2005; but see Dalai et al. 2004; Hennawi et al. 2006; 
Horesh et al. 2005). It is unclear whether these "discrepan- 
cies" are serious. Hence, the systematics in theoretical pre- 
dictions need to be explored, and the effect of baryons is one 
of these. 

Baryons can cool via dissipative processes while DM is col- 
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lisionless, and so their density evolutions may be quite differ- 
ent. Although the baryonic matter is only a modest fraction 
of the total mass (less than or equal to the universal baryon 
fraction of ~ 16%, as revealed by the Wilkinson Microwave 
Anistropy Probe [WMAP]; Spergel et al. 2003, 2006), it can 
nevertheless play a significant role in reshaping the density 
profiles of the DM halos. For example, the baryons that con- 
dense toward the center of DM halos can compress DM by 
adiabatic contraction (Blumenthal et al. 1986; Mo et al. 
1998). This effect was confirmed in simulations with gas 
cooling where the DM density profiles appear steeper than 
those in the pure DM simulations (e.g., Pearce et al. 2000; 
Gnedin et al. 2004). Even for nonradiative gas simulations, 
several studies (Pearce et al. 1994; Navarro et al. 1995; Rasia 
et al. 2004; Jing et al. 2006) suggest that the mass density 
profiles of DM halos become steeper. Some analytic stud- 
ies (e.g., Zhan & Knox 2004) often adopt DM profiles (such 
as the NFW profile) straight from DM-only simulations and 
do not take into account the modification of the DM profiles 
by baryons. In this paper, we explore the effects of baryons 
using high-resolution simulations with nonradiative gas. In 
this case, the interaction between the hot gas and DM is a 
well-defined problem, and hence our study sets a benchmark 
for the effects of baryons on the matter distribution when im- 
portant (but more uncertain) star formation and feedback pro- 
cesses are incorporated. To do this, we perform a set of N- 
body/SPH simulations with different mass and force resolu- 
tions and compare the results with those in a pure DM sim- 
ulation. We will examine effects of resolutions and spurious 
two-body heatings between the gas and DM particles (Stein- 
metz & White 1997). The layout of this paper is as follows. 
In §2 we give an introduction of our simulations and methods; 
in §3 we present our main results. We examine the hypothe- 
sis of energy transfer from DM to gas in §4, and finish with a 
summary and discussion in §5. 

2. THE SIMULATIONS AND METHODS 

We use the massively parallel GADGET2 code (Springel 
et al. 2001; Springel 2005) to simulate structure forma- 
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TABLE 1 
Simulation parameters 



Run Name Num Ngas «jdm "Jgas Softening Length 
/r'Mp r'M Q /T 1 kpc 



PDM 


512 3 - 5.5 x 10 s 


4.5 


Al 


256 3 256 3 3.7 X 10 9 7.4 X 10 s 


4.5 


A2 


256 3 256 3 3.7 X 10 9 7.4 X 10 8 


9.0 


A3 


512 3 256 3 4.6 x 10 8 7.4 x 10 8 


4.5 


A4 


512 3 512 3 4.6 x 10 8 9.2 x 10 7 


4.5 



tion. The code can follow a collisionless fluid with the iV-body 
method, and an ideal gas by means of smoothed particle hy- 
drodynamics (SPH). The GADGET2 implementation of SPH 
conserves energy and entropy (Springel & Hernquist 2002). 
The simulations were performed in a concordance cosmologi- 
cal model with the following parameters (Q m ,Sl\,Sli,,as,n,h) 
= (0.268,0.732,0.044,0.85,1,0.71). Five simulations were 
run with the same initial condition in a cubic box of side- 
length 100 h~ x Mpc (see Table 1 for details). A pure DM 
(PDM) simulation is used as a comparison for the simulations 
that include baryons. In the four simulations with gas, no ra- 
diative cooling is considered and the gas is treated as an ideal 
gas with an adiabatic index 7 = 5/3. 

The highest resolution simulation we performed is A4, 
which includes an equal number (512 3 ) of gas and DM parti- 
cles. The simulation Al has the lowest resolution, but it has 
the smallest softening length relative to the mean interparticle 
separation, about ~ 1 /80 (~ 1 /40 for other simulations). Al 
is used to check the effects of spurious large accelerations in 
close approaches between particles (the real mass distribution 
is inherently smooth and unable to generate large accelera- 
tions). In run A2, we double the softening length from Al 
to check its effects. We also ran a simulation A3, in which 
gas and DM particles have comparable mass to study spurious 
two-body heating effects, which depend on the number of DM 
particles used (Steinmetz & White 1997). Two-body heating 
occurs when DM particles collide with gas particles, lose their 
kinetic energy to the gas components, and, as a result become 
more concentrated. This effect is cumulative in a cosmologi- 
cal simulation and may be particularly serious in simulations 
with low mass resolutions or for halos with a small number 
of particles. Thus it has been suggested that it is wise to use 
more DM particles than gas particles for cosmological sim- 
ulations (Steinmetz & White 1997); the A3 run fulfills this 
requirement. Table 1 lists the number and mass of particles, 
and the softening length used in the simulations. 

The simulations discussed above are all evolved from red- 
shift z = 120 to the present epoch. Because the initial con- 
ditions are identical, one can match each well resolved halo 
on an one-to-one basis in different runs. Thus it is straight- 
forward to select the massive mass halos at z = and com- 
pare their radially averaged density profiles in different sim- 
ulations. The NFW form (Navarro et al. 1997) is used 
to fit the halo density profiles. The NFW profile is given 
by p(x) = Apf(c)/cx/(l + cx) 2 , where A = 200 is the over- 
density, p is the mean universe mass density, x = 7 r 200, c 
is the concentration parameter defined as c = r2oo/r s , and 
f(c) = c 3 /3[ln(l +c)-c/(l + c)]. Here r2oo is the virial ra- 
dius within which the average mass density is 200 times the 
mean background density (not the critical density as used in 
some other studies), and r s is the scale radius. 

The radii of the halos are divided in uniform logarithmic 



steps from 0.001r2oo to 2r2oo- The number of particles and 
the average density in each bin are then computed. Note that 
when fitting the NFW profiles, we use only data between a 
minimum radius, r min , and r 2 oo- For runs Al, A2 and A4, the 
minimum fitting radius r m ; n is 0.02r2oo (corresponding to a 
physical radius from 10.5 to 45 h~ l kpc); for run A2 we use 
'"min = 0.04r2oo to match the larger force softening length. The 
r m ; n value is chosen correspondingly for the PDM simulation. 
The best-fit parameters (including the concentration parame- 
ter c) are found by minimizing 

^2 [log p(n) -log p NFW (n)] 2 . (1) 

i 

Notice that when we fit the DM profiles in simulations with 
gas, we scaled the overdensity A by the DM mass fraction 
(0.84). 

3. THE INFLUENCE OF NONRADIATIVE GAS ON DENSITY 
PROFILES 

We choose a large number of massive halos (about 400) to 
statistically quantify the systematic influence of the hot gas 
on the matter distribution of halos. We only use those halos 
more massive than 10 13 /i _1 M (corresponding to <~ 1.67 x 
10 4 particles in the PDM simulation and 3.11 x 10 4 particles 
in run A4). 

The concentration parameters for these halos in our high- 
est resolution run, A4, versus those in the PDM simulation 
are shown in Figure 1 . The left and right panels show the con- 
centration parameters for the DM and total mass distributions, 
respectively. For the total mass distribution, the concentration 
parameter is increased only slightly (by about 3%), while the 
increase for the DM distribution is larger, at about 10%. The 
effects of baryons (without radiative cooling)-on the density 
profile appear to be modest. 

To examine how the numerical setups affect the DM and 
total mass distributions, in Figure 2 we show the ratios of the 
concentration parameters in nonradiative gas simulations rel- 
ative to the PDM simulation. The median values of this ratio 
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FIG. 1 . — Halo concentration parameters in the highest resolution run, A4, 
compared with those in the PDM simulation. In each panel, the filled circles 
and squares are for the DM and total mass, respectively. The solid lines 
indicate equality, while the dashed lines indicate the concentration parameters 
that are 10% (left) and 3% (right) higher than those in the PDM case. 
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are plotted in five mass bins with equal numbers of halos, to- 
gether with the vertical error bars showing the lower (25%) 
and upper (75%) quartiles in each bin. The ratios of the halo 
concentration parameters in runs Al and A2 with those in the 
PDM run appear to show a systematic trend with mass for 
halos with M < 3 x 10 13 /i _1 M Q . Al and A2 both have 256 3 
DM particles, a factor of 8 lower than A3 and A4. The in- 
ner halo profiles are less well resolved in Al and A2. Thus, 
this poorer resolution results in the mass dependence and the 
higher median values of the ratio of the concentration param- 
eter. The median values for the ratio of the DM in Al and A2 
are comparable with each other. For the Al run, as a small 
softening length (relative to the mean inter-particle separa- 
tion) was used, one can imagine that artificially large accel- 
erations at close encounters may suppress structure formation 
at high redshift in this case, but apparently the effect has not 
accumulated significantly to the present epoch. 

A3 and A4 have the same number of DM particles and 
force resolution; they only differ in the number of gas par- 
ticles (256 3 vs. 512 3 ). The third and fourth columns in Figure 
2 show that A3 and A4 are in good agreement with each other. 
It appears to make little difference for the density profiles of 
halos whether an equal mass is used for DM and gas parti- 
cles in simulations. Furthermore, for both the DM and total 
mass the ratios of their concentration parameters to those in 
the PDM simulation do not have significant dependence on 
the halo mass (especially for A4). The very weak mass de- 
pendence and larger scatter in run A3 may be due to the low 
mass resolution of the gas component. The consistency be- 
tween A3 and A4 shows that the resolutions of A4 are suf- 
ficient for examining the influence of baryons on the density 
profiles for massive halos that are well resolved by a large 
number of particles. 

The results we have shown so far are statistical; below we 
use two individual halos to see the changes in the halo pro- 
files. Figure 3 shows the density profiles, r 2 p(r), for two halos 
with mass of 7.08 x 1O 14 M and 4.16 x IO^/t'Mq. The first 
halo is the most massive one in our simulation. It has a virial 
radius of 2.2/T 1 Mpc and contains 1.3 million DM and 1.2 
million gas particles in run A4, so it is well resolved. Figure 3 
shows that in run A4, the DM and total mass distributions are 
slightly steeper than the DM distribution in the PDM run, be- 
tween 0.02r2oo and 0. 12r2oo, making the concentration param- 
eter larger. As expected, the distribution of gas is substantially 
more extended due to pressure that resists the gravitational 
collapse. This is also reflected in the baryon fraction as a 
function of radius: in the central region, the gas mass fraction 
is only about 10% at R ~ 0.02r2oo (and even smaller at smaller 
radii) while it approaches the universal baryon fraction close 
to the virial radius(see also Ascasibar et al. 2003). For the 
group-sized halo shown in the bottom panel of Figure 3 (other 
low-mass halos show similar behavior), all the simulations 
have roughly the same gas density profiles when R > 0.1r2oo, 
but there appears to be a weak trend for A4 to have a slightly 
lower gas density at intermediate radius (r ~ 0.2r2oo) than Al 
and A3. However, in the inner regions, the gas density distri- 
butions differ markedly: the highest resolution run, A4, has 
the highest gas density, while the lowest resolution run, Al, 
has the lowest values. 

In Figure 4, we plot the temperature and entropy profiles 
for the two halos shown in Figure 3. The trend of the temper- 
ature as a function of radius is roughly the reverse of the trend 
for the density, i.e., A4 has the lowest temperature in the in- 
ner regions. This is because the gas pressure (P oc pT) is well 



constrained by hydrostatic equilibrium (see Fig. 16 in Frenk 
et al. 1999). As one can see in the top panels, the temperature 
profiles exhibit a flat core toward the halo center, and drop 
off sharply toward the outer part of halo. These behaviors 
are similar to those found by previous studies (Loken et al. 
2002; Ascasibar et al. 2003; Rasia et al. 2004) in spite of the 
higher mass and force resolutions in our simulations. Com- 
parisons between different simulations also reveal that as the 
resolution increases, the central temperature decreases while 
the temperature at large radii (r > 0.25r2oo) remains roughly 
the same. Our predicted temperature profiles are somewhat 
different from the observed ones ((e.g., Vikhlinin et al. 2005)) 
which often show a similar behavior in the outer region but a 
decrease in the innermost region. The difference likely arises 
due to radiative cooling, which is not accounted for in our 
nonradiative simulations. 

In the bottom panels of Figure 4, we plot the profiles of en- 
tropy defined as kT /n 2 / 3 , where n is the number density of 
particles. These entropy profiles are consistent with previous 
results (e.g., Frenk et al. 1999; Ascasibar et al. 2003; Kay 
2004; Borgani et al. 2004; Voit et al. 2005). In particular, 
the entropy floor in the central region agrees well with that 
found byAscasibar et al. (2003) who also performed their 
simulations with GADGET2 (see also Voit et al. 2005). In the 
outer region of halos, different SPH simulations yield consis- 
tent results. However, in the inner regions, we find that the 
higher the resolution, the lower the entropy floor. This can 
be understood as follows: at constant pressure, the entropy 
s oc r 5//3 ; in the simulations with lower resolution the temper- 
ature is higher; as a result, the entropy will also be higher in 
low resolution simulations. Note results with the same SPH 
resolution but different DM resolutions (i.e., run Al and A3) 
for the massive halo shown in the bottom left panel, but for 
the low-mass halo, the temperature and entropy in the very 
inner region in run Al are somewhat higher than those in run 
A3. 

Two-body heating is expected to be serious for poorly re- 
solved halos with few particles: gas particles gain energy from 
collisions with DM particles and the gas distribution expands 
substantially. The trend of differences in the gas density pro- 
files in low-mass halos for simulations with different resolu- 
tions is in agreement with this two-body heating effect. Stein- 
metz & White (1997) gave a rough estimate for when this 
effect is important (see their eqs. [5] and [6]). Even for the 
low mass halos shown in Fig. 3, the total number of DM 
particles (AO exceeds 9000 in Al and 76,000 in A4. Taking 
In A = 5 in their equation (6), the two-body heating timescale 
is roughly 200A^/10 4 times the crossing time at half-mass ra- 
dius, which may be too long (see also the left panel in their 
Fig. 2, which shows negligible two-body heating for a galaxy- 
sized halo with 4000 particles). At least for the halos with 
mass larger than 10 13 h~ l Mq in our simulations with highest 
mass resolution (run A4), two-body heating should be negli- 
gible. 

4. ENERGY TRANSFER FROM DM TO GAS PARTICLES IN 
MERGERS 

It is not entirely clear what causes the steepening in the DM 
profiles. One possibility is that during the assembly process, 
shocks produced during mergers between the progenitors can 
convert kinetic energy into thermal energy of gas particles. To 
conserve the total energy, DM particles lose some of their en- 
ergy and sink further toward the center, which steepens the 
DM density profiles. As the merging histories of clusters in 
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FIG. 2. — Mass dependence of the ratios of halo concentration parameters in runs Al -A4 relative to those in the PDM run. The top and bottom panels are 
results for the DM and total mass, respectively. The symbols are the same as Fig.l. The solid horizontal lines indicate equality. The median values and lower 
and upper quartiles in each bin are indicated for five mass bins with equal numbers of galaxies. The median value and the range of mass in each bin are shown as 
horizontal error bars. The number of halos is 399, 400, 395, and 410, for the Al, A2, A3, and A4 runs, respectively. 



cosmological simulations are complicated, they are not ideal 
for isolating the physical mechanism. Instead, we use non- 
radiative mergers between two spherically symmetric clus- 
ters with NFW density profiles to demonstrate the energy ex- 
change between DM and gas particles. 

We utilize the simulations of McCarthy et al. (2006). These 
simulations are gas-only and gas+DM simulations of head-on 
mergers with two equal-mass clusters using GADGET2. In 
the latter case, the gas traces DM initially and the ratio of 
gas to DM mass is set to f4/f2 m « 16%. The virial mass 
of each cluster is 10 15 M Q within r2oo- In these simulations 
each cluster has 50, 000 gas particles and 50, 000 DM particles 
when present. The two clusters are initially well separated, 
with a distance much larger than r2oo- The mergers start with 
the same infall velocity and separation in the gas-only and 
gas+DM simulations, so the initial masses and total energy 
are the same. The force-softening length is set to lOkpc for 
both the gas and DM particles. More details, including the 
setup of initial conditions, can be found in McCarthy et al. 
(2006). 

Figure 5 shows the resulting binned gas density profiles 
from the gas-only and gas + DM simulations after evolution 
for a Hubble time. Both the results of the gas-only and gas + 
DM simulations are similar outside r ~ 150 kpc. But inside 
this radius, the gas density is higher in the gas + dark mat- 
ter simulation, implying that the dark matter must have trans- 
ferred energy to the gas inside this radius. Figure 6 shows 
the evolution of the total energy of the DM from the gas + 
DM simulation. The energy has been normalized by its ini- 
tial value. Note that the total energy is negative (the system 
is gravitationally bound), so if the ratio exceeds unity this im- 
plies that energy has been transferred to the gas. The figure 
clearly demonstrates that at the end of this simulation, approx- 
imately 7% of the DM's total initial energy has been trans- 
ferred to the gas. The same energy transfer procedure should 
have also occurred in the cosmological simulations. 

5. SUMMARY AND DISCUSSION 

In this paper, we performed a set of numerical simula- 
tions with nonradiative gas in the concordance cosmology, 



and compared the results with those of a pure DM simulation. 
The simulations have identical initial conditions but different 
force and mass resolutions. We carefully examined the ef- 
fects of force resolutions and two-body heating. In summary, 
in simulations with hot gas, the DM distribution is more con- 
centrated than those in DM-only simulation. The DM and 
total mass concentration parameters in the nonradiative gas 
simulations are on average about 10% and ~ 3% larger than 
those in the PDM simulation. We used idealized simulations 
of mergers of two clusters to demonstrate that the steepening 
of the DM profile is due to energy transfer from DM to gas 
particles (in shocks), which should also occur in our cosmo- 
logical simulations. 

The influence of baryons on the total mass distribution in 
our simulations is less significant than the results by Rasia et 
al. (2004) who reported a 10% increase in the concentration 
parameter for the total mass distribution. The difference is 
likely due to different SPH treatments used; in GADGET2 the 
entropy is conserved, while in GADGET1 it is not. New re- 
sults obtained using GADGET2 by the same group are consis- 
tent with ours (A. Rasia 2006, private communication). Note 
also that we have many halos (about 400) while they had only 
17 halos, so our statistics are better. Moreover, Ascasibar et 
al. (2003) also showed (their Fig. 1) that gas densities at the 
inner part of halos are higher in a conventional SPH code than 
those obtained in an entropy-conserving SPH code. Overall, 
the influence of hot gas on the mass distribution is quite weak 
in nonradiative gas simulations, especially for the total mass. 
Observational measurement of this effect will be challenging. 
For example, the small change in the total mass concentration 
parameters will not have a substantial impact on the cross sec- 
tion of giant arcs. 

The gas density profiles in our simulations show consider- 
able scatter in their inner parts in different simulations, well 
beyond the softening length (see Fig. 3). We also show the re- 
sults of temperature and entropy profiles to explain the effects 
of the resolution of SPH particles. While two-body heating of 
gas particles by DM particles is qualitatively consistent with 
the trend we see, quantitatively the timescale may be too long 
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for two-body heating to have substantial effects on the gas 
density profiles, in particular for the halos studied in run A4. 
Notice, however, that the uncertainty in gas profile has little 
effect on the DM, and total mass distributions in our results as 
gas only contributes a small (but somewhat uncertain) fraction 
of the total mass. 

The simulations presented above neglected important phys- 
ical processes such as star formation, supernovae feedback, 
and heat conduction. For the group- and cluster-scale halos 
that we analyzed, the processes mentioned above are expected 
to be less important than in galaxies, nevertheless, their proper 
treatment is important in comparing data and theoretical pre- 
dictions. Our study also ignores the differences in dynamical 
frictions between galaxies and DM halos. Because galaxies 
are located at the centers of DM halos, they can survive the 
tidal disruptions longer than the DM halos (e.g., Gao et al. 
2004). As the dense cores of galaxies sink toward the cen- 
ters, their kinetic energy can be transferred to and heat up the 
DM particles to make the DM distribution flatter (El-Zant et 
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FIG. 3 . — Density profiles for the most massive halo (top) and a group-sized halo (bottom) in our simulations. The rilled dots indicate the density profile for the 
PDM run, while other lines indicate results for the total, DM and gas density profiles in runs Al, A3, and A4 (labeled inside the panels). The three dotted vertical 
lines indicate the softening length, 0.02r2oo and rioo- The data points between the latter two are used in the profile fitting in Al, A3, and A4. 
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FIG. 4. — Temperature (top) and entropy profiles (bottom) for the two halos shown in Fig. 3. The line symbols are the same as in Fig. 3. 
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FIG. 5. — Gas density profiles in the gas-only (solid line) and in the gas + DM (dashed line) simulations of head-on collisions between two spherical NFW 
clusters (see §4 for more details). 



T | I | I | I | I | I [■ 



1.1 - 




2 4 6 8 10 12 

time (Gyr) 



FIG. 6. — Ratio of the total energy of the DM at time t to the initial total energy in the gas + DM simulation. The dotted line indicates energy conservation. At 
t > 5 Gyr, the dark matter has lost energy to the gas, as implied by a ratio larger than unity. 



